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Particles interacting with short-ranged potentials have attracted increasing interest, partly for 
their ability to model mesoscale systems such as colloids interacting via DNA or depletion. We 
consider the free energy landscape of such systems as the range of the potential goes to zero. In this 
limit, the landscape is entirely defined by geometrical manifolds, plus a single control parameter. 
These manifolds are fundamental objects that do not depend on the details of the interaction 
potential, and provide the starting point from which any quantity characterizing the system - 
equilibrium or non-equilibrium - can be computed for arbitrary potentials. To consider dynamical 
quantities we compute the asymptotic limit of the Fokker-Planck equation, and show that it becomes 
restricted to the low-dimensional manifolds connected by "sticky" boundary conditions. To illustrate 
our theory, we compute the low-dimensional manifolds for n < 8 identical particles, providing a 
complete description of the lowest-energy parts of the landscape including floppy modes with up to 
2 internal degrees of freedom. The results can be directly tested on colloidal clusters. This limit 
is a novel approach for understanding energy landscapes, and our hope is that it can also provide 
insight into finite-range potentials. 



The dynamics on free energy landscapes is a ubiqui- 
tous paradigm for characterizing molecular and meso- 
scopic systems, from atomic clusters, to protein folding, 
to colloidal clusters. [IHl]- The predominant strategy 
for understanding the dynamics on an energy landscape 
has focused on the stationary points of the energy, the 
local minima and the transition states, and seeks the dy- 
namical paths which connect these to each other, while 
more recent models generalize to metastable states con- 
nected by paths as a Markov State Model [5]. These 
techniques have proven to be extremely powerful, giv- 
ing innumerable insights into the behaviour of complex 
systems j6HT3], On the other hand, a major issue has 
been the difficulty of finding the transition paths, con- 
necting local minima or metastable states to each other, 
especially given a complex energy landscape in a high- 
dimensional space. A variety of creative methods have 
been developed in recent years for efficiently finding tran- 
sition paths p!4H24] but for a given system, there is no 
guarantee that all such paths have been found. 

Here we present a different point of view for under- 
standing an energy landscape, that occurs when the 
range over which particles interact is much smaller than 
their size. Such is the case in certain mesoscale sys- 
tems, for example, for Ceo molecules 25l[26], or for col- 
loids interacting via depletion [4 or coated with poly- 
mers or complementary DNA strands [27H3Q]. We will 
show that in this limit, the free energy landscape is de- 
scribed entirely by geometry, plus a single control pa- 
rameter n that is a function of the temperature, depth, 
and curvature of the original potential. This limit is 
related to the sticky sphere limit of a square-well po- 
tential [31], which has been used to investigate thermo- 
dynamic properties of hard sticky spheres j32H34j . The 
landscape can be thought of as a polygon living in a high- 



dimensional space, whose corners (0-dimensional mani- 
folds) are joined to each other by lines (1-dimensional 
manifolds), that in turn form the boundaries of faces (2- 
dimensional manifolds), and so on. These manifolds are 
fixed functions of the particles' geometries, independent 
of the details of the original interaction potential from 
which the limit was taken. 

Once the geometrical manifolds comprising the land- 
scape are computed, non-equilibrium quantities charac- 
terizing the dynamics can be calculated by solving the 
Fokker-Planck equation or its adjoint on these manifolds. 
We show that in the short-ranged limit these equations 
acquire an effective boundary condition at the boundary 
of every p-dimensional manifold in the polygon. This 
makes the kinetics computationally tractable, since the 
stiff modes of a narrow potential become a set of bound- 
ary conditions. 

The geometrical nature of the energy landscape does 
not mitigate its high dimensionality, but at low tem- 
peratures (high n) both the free energies and the kinet- 
ics are dominated by the lowest-dimensional manifolds. 
This means that the description of the free energy land- 
scape and kinetics for short-ranged potentials reduces to 
a problem in discrete and computational geometry - to 
find all of the low dimensional manifolds for a given set 
of interacting particles. 

As an illustration of the framework, we characterize the 
geometrical landscape for n < 8 particles with identical 
potentials, and demonstrate how these solutions lead to 
a complete description of the energy landscape and the 
kinetics of this system. This solution describes both the 
geometrical limit of small atomic clusters as well as a di- 
rect prediction for colloidal clusters interacting with de- 
pletion forces [H ini EH |36] , where the predictions could 
be tested experimentally. The solution also provides a 
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framework for understanding and extending simulations 
on clusters with finite range potentials [1 . Our calcula- 
tion of the energy landscape builds on the enumeration of 
all finite sphere packings of n particles with at least 3n — 6 
contacts [371 EH]- With these as the starting point, we 
compute every 1- and 2-dimensional manifold of motions 
that contains a 0- dimensional manifold at its boundary, 
from which we can extract statistical quantities such as 
the relative entropies of the different types of motions. 
Then, we solve Fokker Planck equations on these mani- 
folds to obtain transition rates between the lowest energy 
states, the 0-dimensional manifolds. 



THE GEOMETRICAL LANDSCAPE 

We begin by showing how the geometrical free energy 
landscape arises as a distinguished limit of particles in- 
teracting with arbitrary potentials of finite range. We 
consider a point x G in configuration space, with 
the potential energy given as a sum of potentials concen- 
trated near different geometrical boundaries as 



U'{x) = ^CMyk{x)/e). 



(1) 



The functions j/fe : R" ^- R, (fc = 1, . . . , 

kmax) represent 

the geometrical boundaries via their level sets {yk{x) = 
0}, and U{y) : R ^ R is the potential energy near 
each boundary. For concreteness, let us suppose this is a 



model for n spheres with centers at Xi G 



1 . . . n. 



so the configuration is x = (xn, X12, X13, X21, . . . , Xns) G 
R^"^, and we take i/k to be the excess bond distance, as 
i/k = ^/\xi^k)'^^^^jik)\^ — where d is the particles' di- 
ameter. Then U {y) is the pairwise interaction potential 
that we assume has minimum Uq dit y = and is negli- 
gible beyond some cutoff Tc- For ease of exposition, the 
interaction potential is assumed to be identical between 
each pair of particles, but this is not a necessary restric- 
tion for the geometrical landscape to apply. For the total 
potential t/^, the parameter e characterizes the range of 
the potential, while Ce is proportional to the depth. The 
geometrical limit requires taking Ce 00 as e in a 
manner that we specify momentarily. 

We consider particles that evolve according to the over- 
damped Langevin dynamics [39| at temperature T: 



dx 
'dt 



1 

7 



VU'{x)^V2Dr]{t), 



(2) 



where 7 is the friction coefficient, D = is the 

diffusion coefficient, /3 = {kBT)~^^ /c^ is the Boltzmann 
constant, and 77 (t) is a 3n-dimensional white noise. The 
equilibrium probability for this system is the Gibbs dis- 
tribution [40 : 



dp'{x) = {Z')-^e-^^^^^Ux. 



(3) 



where is the normalizing constant. 

The geometrical free energy landscape occurs when 
the range e ^ 0. The relationship between the depth 
and the range is critical to obtaining an interesting 
limit. If only the range goes to zero, then particles 
are bonded for shorter and shorter times so that in 
the limit they simply behave like hard spheres. On 
the other hand, if the depth goes to —00 too quickly, 
then the particles simply stick together and only un- 
bind on exponentially long timescales [41 . The inter- 
esting limit occurs when particles stick to each other 
but unbind on accessible timescales; for this we must 
choose Ce so that the Boltzmann factor for two particles 
to be bonded to approaches a finite, non-zero constant: 
K, = lime^o ^ Jq^"" e~^^^^^^/^^(ix, where we define Boltz- 
mann factors non-dimensionally by scaling by the diam- 
eter d. Evaluating the integral using Laplace's method 
then implies 

n = lim — (4) 

d^/c(3CeU^\0) 

We call the constant the sticky parameter. Note that hc 
is a function of both the potential and the temperature. 
The constant c = 2/7r for hard-spheres. 

In the geometrical limit, the probability measure 
becomes concentrated at the exact locations in config- 
uration space where a bond forms, i.e. on the level 
sets {yk{x) = 0} and all possible multi-way intersec- 
tions. Thus the limiting probability distribution will be a 
weighted sum of delta functions, each defined on a mani- 
fold corresponding to a different set of bond constraints. 
The weight of the each delta function depends on the 
number of bonds and a geometrical factor associated with 
the entropy of the configuration. This gives a geometri- 
cal picture of the energy landscape: When hz is large, the 
occupation probabilities will be dominated by configura- 
tions where the number of bonds m is large. For identical 
particles, with n < 9, the maximum number of contacts 
is m = 3n — 6 [371 EH]' these are rigid structures that 
have no internal degrees of freedom, so they correspond 
to 0-dimensional manifolds, or "points". The next low- 
est configurations in potential energy are manifolds with 
m = 3n — 7, which are are obtained from rigid structures 
by breaking one bond. These have one internal degree of 
freedom so are 1-dimensional manifolds, or "lines". The 
lines form the boundaries of 2-dimensional manifolds, or 
"faces" , when another bond is broken, and continuing up 
in dimensionality we obtain the entire energy landscape 
as the union of manifolds of different dimensions. 

Figure [T] shows a schematic contrasting this limiting 
case with the traditional picture of an energy landscape. 
The traditional picture is of an undulating surface, with 
local minima connected through saddle points, whose 
heights provide an activation barrier that determines the 
transition rates between basins. In contrast, in the ge- 
ometrical limit, the local minima are infinitely narrow 
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and deep, separated by long, relatively fiat spaces in be- 
tween - the landscape can be thought of as a golf course 
punctuated with deep trenches very deep holes. Kinetics 
on this landscape are determined partly by an activation 
barrier - the time it takes to climb out of the hole - and 
partly by diffusion. 

The figure also shows an explicit 1-dimensional mani- 
fold taken from the landscape for n = 6, an example we 
will return to throughout the text. There are two ground 
states, the polytetrahedron and the octahedron [38 , and 
the manifold is the set of deformations corresponding to 
the transition path between these when a single bond is 
broken. 

To explicitly calculate the equilibrium probabilities of 
the different states in the geometrical landscape we con- 
sider a configuration with m constraints or equivalently 
p = 3n — 6 — m bonds broken. The constraints are writ- 
ten as an ordered multiindex a = (ai, 0^2, . . . am) so the 
manifold of configurations satisfying such constraints is 

^a={x : Va^ {x) = VaA^) = . . . = {x) = 0, 

y^{x) >0,/3^ai}. (5) 

We write a = to mean the region where no constraints 
are active, and let Q = Ua^a be the full space of ac- 
cessible configurations. The limiting partition function 
associated with these constraints is 

= lim / e-^^^^'^Ux, (6) 

e^O d^n V ; 

where is the neighbourhood surrounding the manifold 
where the potential associated with the constraints is 
active: 

ni^ = {x :0< ya^{x),...,ya^{x) < ere, 

y^{x) > er^P ^ ai}. (7) 

This splits configuration space near each manifold into 
two parts - the fast variables y^y. changing rapidly along 
directions associated with the constraints (sometimes 
called the vibrational modes), and the slow variables yj^ 
that are the unconstrained configuration. 

To compute the integral in ^ we need a parameteriza- 
tion of the manifolds associated with the constraints. It 
is convenient to parameterize the fast directions by the 
constraint variables themselves, x Vai^ • • - yam- We 
choose the additional 3n — m variables y G R^^-^ so 
that Vy-Vyai = on I^q,, i.e. the variables y^ y^^ are or- 
thogonal on the manifold. As discussed in the Appendix, 
it is possible to find such a parameterization locally as 
long as the manifold associated with the constraint vari- 
ables is regular - i.e., the coordinate transformation for 
the constraint variables must be smooth and invert ible. 
This happens when the Hessian of the potential energy 
(say at e = 1) 

m 

H^{x) = VVU'=Hx) = f/"(0)^Vy„.(Vy„J^ (8) 



has m non-zero eigenvalues. 
We can now evaluate 



e e 



-/5a EI^i UiYi) 



(9) 



where Y = (Yi, . . . , Y^) with Yi = y^./e, gij is the metric 
tensor associated with the transformation (^q, . ^y) ^ x 
with components gij = JkiJkj = J^J, where Jij = 
and \gij\ is its determinant. Let separate the metric ten- 
sor into blocks as 
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9ab 9av 
9ub 9uv 



(10) 



The first set of indices run from 1 < a, 6 < m and 
describe the fast variables, perpendicular to the mani- 
fold, while the second set run from m -\- 1 < u^v < N 
and describe the slow variables, along the constraint 
manifold. Let \gcd\ be the determinant of a particular 
block. It follows from the definition of the metric that 
|^a6|y=o = ni^i-^r^' where the A^'s are the non-zero 
eigenvalues of Ha{x)/U"{0), and the condition of or- 
thogonality gives \gav\Y^o = \9ub\Y=o = 0. Evaluating 
the integral in ([9| over the fast variables using Laplace's 
method, and letting e ^ 0, shows the limit is 



where 



-1/2 



(11) 



(12) 



is a geometrical factor (representing the "vibrational" de- 
grees of freedom) that depends only on the level sets of 
the constraints, and g^ = guv\Y=o is the metric on man- 
ifold Qa inherited from the ambient space by restriction. 



The integral (11) is simply the volume integral of ha{y) 
over 

The manifold contains 6 degrees of freedom repre- 
senting translation and rotation of the cluster, and the 
partition function integral can be further simplified by 
integrating over the subspace spanned by these motions. 
This introduces a factor I{x) in the integral, the square 
root of the moment of inertia tensor [40 . If we let n2 be 
the quotient space formed by identifying points x ^ z if 
X can be mapped to z by a combination of translations 
and rotations, i.e. — ^a/SE{3), where ^£^(3) is the 
special Euclidean group, then we can write 



/ ha{x)I{x)^/\g^\d2 



(13) 



where g^ is the metric on i we have dropped con- 
stants (such as free volume) that are the same for all 
configurations. For convenience later on, we define (a 
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as the part of the partition function that is independent 
of hv. The Appendix has a detailed discussion for how 
to construct which is critical for using the formalism 
developed here for practical calculations. Fig. [l] (bot- 
tom) is an example of such a quotient manifold, where 
each point on the manifold represents the 6-dimensional 
space of clusters in configuration space that are related 
by rotations or translations. 

The particles in Figure [l] are different colours to iden- 
tify the different transitions that occur when moving 
around configuration space. However, when (as imagined 
here) the particles are identical, permuting the colors of 
any particular structure yields a geometrically isomor- 
phic structure on a separate part of the quotient space. 
The free energy must account for the number of dis- 
tinct manifolds that are geometrically isomorphic to • 
When p = this is = CoN\/a^ where a is the symme- 
try number, i.e. the number of particle permutations that 
are equivalent to a rotation, and Co = 2 if the structure 
is chiral and Co = 1 otherwise [42 . For p > 0, we count 
the multiplicities by counting how many times a mode oc- 
curs from the perspective of each 0-dimensional "corner" 
of the mode, and dividing by the total number of corners. 
[43] For example. Fig. [l] has corners from two different 
ground states, the polytetrahedron and the octahedron, 
which occur with multiplicities ni,n2 respectively. For 
each polytetrahedron, there is ui = 1 line coming out of 
it that is isomorphic to this transition, and for the octa- 
hedron there are = 12 distinct lines. (The numbers 
z^i, 1^2 are indicated on the arrows connecting red circles 
to blue circles in Figure |2j where the transition under 
consideration is mode 7.) Therefore the total multiplic- 
ity of the line is {niUi + 712^^2 )/2. Consider a transition 
connecting a polytetrahedron to a distinct copy of itself, 
say mode 5. Here there are ui = 4 such lines connected 
to each polytetrahedron, so the multiplicity of the line is 
niz^i/2. In general, each p-dimensional manifold con- 
tains a total of Tic corners from N < Uc nonisomorphic 
ground states, each ground state having multiplicity n^, 
and such that each single ground state is connected to Vi 
distinct manifolds isomorphic to so the multiplicity 
IS Ua = l^i^iniUi/ric. 

Putting this all together, the total partition function 
of all structures isomorphic to a given constraint mani- 
fold and the free energv of these structures 
are = — /c^T log (ric^Zc^). We can separate this free 
energy into = — mkBT log n — kBTlog{naCa)^ using 
the definition of (a in Equation (13). The first term 
{—mkBTiogn) depends on the temperature, bond en- 
ergy, and width of the potential, whereas the second term 
with Sex = —^B log{naCa) IS entirely geometrical, and es- 
sentially is the entropy of structures corresponding to 
the constraint set l^J. As an example, we have plot- 
ted Fa/kBT along the polytetrahedral-octahedral tran- 
sition in Figure [1] (bottom). This varies smoothly along 
the manifold as /(x), h{x) vary. The endpoints of the 



manifold are the ground states, where the free energy 
changes discontinuously because the number of bonds has 
changed - this causes a jump in the energetic part via m, 
and the entropic part via h^- 

With these results in hand, we can now compare the to- 
tal entropies of floppy manifolds of different dimensions, 
to understand the temperature range in which the differ- 
ent manifolds occur. Let the total geometrical partition 
function of manifolds of dimension p be 



E 

dim(^7g)=p 



Z = ^f^^''-^-PZp. (14) 



Here Zp is independent of the temperature and poten- 
tial, while Z combines everything to obtain the entire 
landscape. Note that lower dimensional manifolds have 
more bonds and thus are favored in the partition func- 
tion when the temperature (or equivalently hz) is small. 
As temperature increases, k shrinks and higher dimen- 
sional manifolds become more highly populated. Even- 
tually clusters will fall apart into single particles. The 
temperature dependence of how clusters fall apart is en- 
coded in the relative sizes of the ZpS. The temperature 
where the landscape transitions from having more p di- 
mensional structures than p-\-l dimensional structures is 
found by solving hcZp = Zp_^i for which gives roughly 
{kBTp)~^ ^ \n Zp-^i/Zp -\- const. 



Free energy landscape for identical particles 

To illustrate our asymptotic calculations with a con- 
crete example we have computed the geometrical man- 
ifolds up to dimension p = 2 for n = 5,6,7,8 identical 
particles with diameter d = 1. To do this, we begin the 
set of clusters with > 3n — 6 bonds derived in Arkus 
et. al. [38] For every rigid cluster we break each sin- 
gle bond in turn and move along the internal degree of 
freedom until we form another bond. This is the set 
of one-dimensional manifolds. For the two dimensional 
manifolds, we break each pair of bonds from the rigid 
clusters in turn, and move along the internal degrees of 
freedom to compute the boundaries, corners, and inte- 
rior of the two-dimensional manifolds (see Appendix for 
details.) This algorithm ensures we have every floppy 
manifold that can eventually access one of the rigid clus- 
ters in our list only by forming bonds, but never breaking 
them. Our analysis makes three assumptions that we be- 
lieve to be true, but await rigorous proof: First, we are 
assuming that the list of clusters in Arkus et. al. |38] 
is the complete set of rigid (0-dimensional) clusters; this 
is true as long as there are no rigid clusters with 3n — 7 
bonds or less, a condition that was not checked [44J. Sec- 
ondly, in the calculations of the entropy of the two dimen- 
sional floppy manifolds we assume that the manifolds are 
topologically equivalent to a disk. The fact that our pa- 
rameterization algorithm works is evidence for this claim. 
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TABLE I. Geometric partition functions Zp, and numbers of different modes, for the set of p = 0, 1, 2-dimensional manifolds 
as the number of spheres n varies. These are geometrical quantities that do not depend on the temperature or potential. The 
total partition function of the p-dimensional manifolds, which includes both of these dependencies via is Zp. 



n 


#0-D 


# 1-D 


#2-D 


Zo 


Zi 


Z2 


Zi/Z^ 


Z2/Zi 


5 


1 


2 


4 


10.7 


73.8 


545 


6.9 


7.4 


6 


2 


5 


13 


36.1 


256 


1140 


7.1 


4.5 


7 


5 


16 


51 


1.1x10^ 


8.5x10^ 


39 X 10^ 


7.6 


4.6 


8 


13 


75 


281 


49x10^ 


396x10^ 


1.87x10^ 


8.1 


4.7 



though we have not proved this rigorously. Thirdly, we 
assume that all floppy manifolds can eventually access a 
rigid mode and are not, for example, circles. We men- 
tion these caveats because although we are confident that 
they do not apply in the low n examples described here, 
it is possible that potentially significant exceptions arise 
at higher n. 

The landscape for n = 6 is shown [45] in Figure |2j 
There are two ground states, denoted by the red circles, 
each with 3n — 6 = 12 contacts; the area of the red cir- 
cles are proportional to the probability of each state, with 
the polytetrahedral ground state ^ 25 times more likely 
than the octahedral. The light blue circles denote the 5 
topologically unique structures that are missing a single 
bond in the ground states - such structures correspond 
to a one dimensional manifold, with continuous deforma- 
tion along the direction of the missing bonds. The yellow 
circles denote the 13 unique structures that are missing 
two bonds from the ground state. Again, the area of the 
circles is proportional to the occurrence probability of 
these modes. These structures correspond to two dimen- 
sional manifolds, with continuous deformations allowed 
along both of the directions between the missing bonds. 
The connections between the different modes are denoted 
by arrows on this figure, with structures missing (say) 2 
bonds generally arising from breaking a single bond from 
several different pathways. 

Each element of a mode with two bonds broken can be 
mapped to a polygon in R^, and these parameterizations 
are also shown in Figure |2j We have chosen one param- 
eterization (mode 18), to illustrate in detail in Figure [s] 
The interior of the manifold represents structures with 
10 bonds, with each point representing a different set 
of coordinates for the particles. The edges correspond to 
structures with 11 bonds, while the corners are structures 
with a full 3n — 6 = 12 bonds. Of the four corners, three 
correspond to polytetrahedra with different permutations 
of the particles, and the final one is an octahedron. The 
one-dimensional edges connecting the corners are possi- 
ble transition paths that can be followed by breaking only 
one bond. 

In general the number of corners varies among modes, 
with no a priori way to determine this without solving the 
full geometry problem. It ranges from 3-6 for n = 6, 7, 
and from 3-7 for n = 8. Many 2-dimensional modes con- 
tain several permutations of a given 0-dimensional mode 



as a corner. 

Table 1 summarizes the partition function data. The 
number of different modes increases combinatorially with 
n, as do the geometrical partition functions Z^. Strik- 
ingly, the ratios Zi/Zq, Z2/Z1 remain virtually con- 
stant as n increases. This implies that the temper- 
ature dependence of the landscape is independent of 
n. Indeed, Figure [4] shows the relative probabilities of 
0,1, 2-dimensional modes as the temperature varies, for 
n = 6, 7, 8, with the yield of a p-dimensional mode given 
by Hp = ^~'^Z^I(j^Z{^ + kZ\ + Z2). As an illustration, 
we have chosen V^lks = -4, {U'\0)/kB){2/7r) = 15, so 
that i^{T) = e^l^ I^XhjT. Because the ratios Zpj^i/Zp 
are nearly the same, these graphs are essentially indis- 
tinguishable for different numbers of particles. Moreover 
the critical temperature for transitioning from mostly 
dimensional structures to 1 dimensional structures 
(^ 1/lnZi/Zo) is quite close to that for transitioning 
from 1 dimensional structures to 2 dimensional struc- 
tures (~ 1/ In Z2/Z1). If this remains true as p increases, 
it would imply that clusters melt explosively at some 
critical temperature, rather than incrementally: clusters 
would occupy either mostly the 0-dimensional modes, or 
a gaseous, no- or few-bond-state, but not the chain-like 
floppy configurations in between. 



KINETICS ON THE GEOMETRICAL 
LANDSCAPE 

We now consider kinetics on the geometrical land- 
scape. The concentration of equilibrium probabili- 
ties on manifolds with varying dimensions also applies 
to non-equilibrium quantities, such as transition rates, 
first-passage times, the evolution of probability density, 
etc. Such quantities can be computed from the time- 
dependent transition probabilities, which for dynamics 
given by ([2| are obtained from the corresponding for- 
ward or backward Fokker-Planck (FP) equations [46] , 
We show that in the geometrical limit, the FP equation 
asymptotically approaches a hierarchy of FP equations, 
one on each manifold of each dimension. These equations 
are coupled to each other, with equations on manifolds 
with dimension p serving as boundary conditions for the 
equations on manifolds with dimension p + 1. 

The idea behind the derivation is quite natural: if a 
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potential is deep and narrow, then it equilibrates much 
more rapidly in the directions along the bonds than along 
a cluster's internal degrees of freedom. Therefore the 
probability density near a cluster with p bonds broken is 
approximately p{y,t)e~^^i^ where y parameter- 

izes the internal degrees of freedom of the cluster. The 
"constant" p{y, t) evolves slowly in the transverse direc- 
tions according to the Fokker-Planck dynamics (which 
includes an "effective" potential that arises from the cur- 
vature of the manifold), and it also changes due to the 
flux of probability out of the p + 1-dimensional manifolds 
for which it forms part of their boundary. This gives a 
hierarchy of coupled FP equations. 

To see how this comes about in detail, we examine so- 
lutions to the Fokker-Plank equation for the evolution of 
the probability density p^{x,t). Given a parameteriza- 
tion of configuration space M.^^ with metric tensor ^, the 
non-dimensional Fokker Plank equation corresponding to 
^ is 



dtP^ = div {p^ grad + grad p^) 



(15) 



with boundary conditions at each level set {yi{x) = 0} 

{p' grad U'^ grad p')-h^'^ = {p'g'^djU' + g'^djp) hf = 

where n^*^ is the outward normal to the boundary. We 
have non-dimensionalized lengths by (i, times by (P /D^ 
and energy by ksT. Away from all boundaries there is 
no force, so in VLq the limiting probability p evolves only 
by diffusion as 



dtP ■ 



div jo, 



with 



Jo 



grad p. 



(16) 



Now consider the evolution near a manifold l^c^, a 
"wall" . The dynamics in the directions orthogonal to the 
wall, where bond distances are changing, are much faster 
than those along it, so near the wall the probability den- 
sity will rapidly approach a multiple of the equilibrium 
probability. Parameterizing the region near the wall 
as {y^i , y) as in the previous section, we obtain 

p'{yo^^ ,y,t)= p(0, y, t)e- ^^=i ^^^^"^)+epie- ^^=1 

(17) 

where epi{y^yai^t) is the correction to the leading order 
formula. This satisfies the matching condition that p^ 
be asymptotically continuous. This ansatz can also be 



derived from a consistent asymptotic expansion of (15) 
after the change of variables Yq,. = ^a^/e. 



Substituting (17) into the Fokker-Planck equation (15) 
gives 



dt 



(e-^^p^ = div (e-^^ gradp^ + 0(e). (18) 



We would like to to integrate out the fast variables so 
we need to separate these from the slow ones. This is 



most conveniently done using the metric g with block 
decomposition ( p!Q|). Ma king the change of variables Yi = 
^Q,./e shows that (A.36) can be written as 



epi) 



-"-g^^d^ip- 



5„(v^e-f^"5™(9„(p + epi 



)) 



gle-^'-g^^'dhip + ep,: 



(19) 



where we abbreviate Ua = Xli^i CeU{Yi). 

We now integrate ([l9| over the fast variables 
e^dYidY2 . . . dY^ and keep the leading-order parts. The 
terms oc ^ vanish in the limit because g^^\Y=o = 
g^^\Y=o = 0- The 0(1) terms require evaluating an inte- 
gral similar to (|9| which converts ^/\g\e~^'^ to the factor 
\/\guv\f^^ha at each point on the manifold. 

The first term is the most interesting. This is the 
divergence in the fast variables, and although p does 
not depend on these, the unknown pi might and con- 
tributes at 0(e~^). However, we can use the divergence 
theorem to replace this term with an integral of the 
flux through the p-dimensional fast-variable boundary, 
at {x : Yi{x) = Tc}. We then introduce a second match- 
ing condition which requires the flux to be continuous, 
so we can replace it with the sum of fluxes from each 
p + 1-dimensional manifold which has Qa as part of 
its boundary, [47] and evaluate the limit of the boundary 
integral for these matched fluxes. 

Finally, we integrate over the space of rotations and 
translations of each point, assuming p(x, 0) is constant on 
orbits. The divergence in these directions will disappear 
by Stokes' theorem, and the remaining directions provide 
dynamics on the quotient space. Combining with the 
previous calculations, yields 



1 



(^V\9uv\l^ag''''dyp^ 



divc, {Ka grade, p) 



where the fluxes have leading order part 

ja = -i^a grade, P onn2- 



(20) 



(21) 



Here hz^ = i^^h^I combines the sticky parameter and 
the geometric factor at the wall, and the metric tensor 
g is the quotient metric obtained from the metric ga on 
which in turn is inherited from the original metric g 
in the ambient space by restriction: g^ = g\Yi=...=Ym,^o- 
The sum is over (3 such that ft^ Part of the bound- 
ary of the p + 1-dimensional manifold , and where 
n^" is an outward normal vector to at Q^- We de- 
fine diYa 7 grade to be the differential operators on the 
quotient manifold ft^- 
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Equation (20) has a more intuitive interpretation as 
the evolution of the probabihty along the wall It is clear 
from the derivation that the total probability density 
(with respect to the wall coordinates) of being on the 
wall is Pa = i^aP- This satisfies 



dtPa 
1 



divc, -Pa grade log i^a + grade, Pcx 



J/3 • n' 



effective force 



diffusion 



The dynamics along the wall is therefore a combination 
of diffusion, plus drift due to an "effective" potential 
— log A^Q,, plus a flux in and out of the wall. 

The effective potential is entropic in nature and comes 
from the changing wall curvature, which makes the po- 
tential look wider in some places than others. A parti- 
cle will spend more time in the wide places than in the 
narrow ones, and since it reaches equilibrium much more 
quickly in the transverse directions than in the along- wall 
directions, it looks like there is an effective force pushing 
it to the wider areas. This is the same equation one ob- 
tains by letting the depth of the potential become infinite 
without changing the width but with the addition 
of the flux in and out of the wall. 

This flux term is illustrated schematically in Figure |5] 
for a simple case where the configuration space (x^y) G 
has a single constraint, yi{x, y) = x - a. "wall" at the 
horizontal axis. The probability integrated over a small 
box (red) with length Ax near the wall changes in a time 
increment At due to two processes: probability fluxing 
along the wall in the x-coordinate, which contributes a 
change of At{px{x-\- Ax^t) —px{x^t)) J e~^^y^dy^ and the 
flux of probability from the wall to the interior, which 
contributes a change of AtAxp^(x, 0, t). Equating with 
Atpt{x,t) J e-^^y^y gives the effective boundary condi- 
tion. 

To summarize, the limiting FP equation and boundary 
conditions are (16), plus (20) on every ft^- Substituting 
for the time derivatives shows that the boundary con- 
ditions are second-order, and this is why conditions are 
needed on every boundary and not only those with co- 
dimension 1. We call this set of equations the "sticky" 
equations because it is the Fokker-Planck equation for 
Sticky Brownian Motion [48^, a stochastic process that 
has a probability atom on the boundary of its domain. 



TRANSITION RATES FOR STICKY BROWNIAN 
MOTION 

To transition between the rigid configurations in the 
geometrical limit, a cluster must break one (or more) 



bonds, then diffuse across the line segment (or face, etc) 
until it hits the other endpoint. The time it takes to do 
this depends on the length of the line and on how the 
entropy of the configuration varies along the line, and 
we can find this by solving an equation on the full line 
segment (face, etc). In our asymptotic limit there are no 
meaningful transition "states" - rather, the entire line 
segment can be thought of as a transition state. Once the 
energetic barrier of breaking a bond has been overcome, 
transitions are dominated by diffusion. 

We now consider how to compute transition rates using 
the sticky equations (20), supposing we have a stochas- 
flux to/from wall tic proccss X (t) whosc probability evolution is well- 
(2i§))proximated by these. [49 Consider the transition rate 
between sets A C fl^ B C fl^ and for simplicity let us 
focus only on the case where these are both (disjoint) 
subsets of the 0-dimensional manifolds. For example, we 
may be interested in the transition rate between the oc- 
tahedron and the polytetrahedron (introduced in Figure 
[T]), in which case A would contain all points in the quo- 
tient space representing the octahedron, and B all those 
representing the polytetrahedron. These rates can be 
computed using Transition Path Theory, which provides 
a mathematical framework for computing transition rates 
directly from the Fokker-Planck equations. We will sim- 
ply state the facts that are relevant to our example and 
refer the reader to other resources for more details [5Ql - 
52 ). 

The committor function q{x) is the probability, start- 
ing from point x, of reaching set B before set A. This 
solves the stationary backward Fokker-Planck equations 
with boundary conditions q{A) = 0, q{B) — 1, plus any 
other boundary conditions remaining from the equations. 
As shown in the Appendix, the forward sticky equations 
(20) are self-adjoint (with respect to the invariant mea- 
sure) so these are also the backward sticky equations. 

A reactive trajectory is a segment of the path X(t) that 
hits B before A going forward in time, and A before B go- 
ing backward in time. The probability current of reactive 
trajectories is a vector field that, when integrated over a 
surface element, gives the net flux of reactive trajectories 
through it. Because our process is time-reversible, this 
current is [51J 



J(x) = p(x)V(7(x), (23) 

where p{x) is the equilibrium probability measure. From 
(13) we find this is 

=Z-i^K„(a;)5„(x), (24) 

a. 

where ^q;(x) is the singular measure on I^q,, i.e. it satis- 
fies J^f{x)Sa{x) = f{x). Finally, the transition rate 
k"^^ is calculated by integrating this flux over a surface 
S dividing the two states, giving 



L 



= I J-ndS, 



(25) 
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where n is the normal pointing from the side containing 
A into the side containing B. 

Computing this exactly requires solving the backward 
equations on a high-dimensional space, and integrating 
over a high-dimensional surface - a computationally in- 
feasible proposition. However, when the sticky parameter 
K is large, most of the probability is concentrated on the 
lowest- dimensional manifolds so we expect these to con- 
tribute the most to the rates. Therefore, let us expand 
the equations in powers of f<i~^. We suppose all variables 
have an asymptotic expansion as 



(26) 



and similarly for ja^ p, J, etc. Expanding p shows that 
to first order it is a measure on points, to second order 
it is a measure on points and lines, etc. Measures on 
points will not contribute to the rate because the dividing 
surface S can be chosen to avoid them, so k^^ = and 
the leading-order part of the rate is 0{i<i~^)^ computed 
from piVqo- 

Expanding the backward sticky equations in powers of 
gives a set of equations for qq: 

dWa hal grade <7o = on > 0), (27) 

with boundary conditions qo{A) = 0, qo{B) = 1, and 
X]/3^ck(^/3^ grad^ qo)-n^^ = on all other 0-dimensional 
manifolds. 

To solve these equations we can first find the solution 
on the 0- and 1-dimensional manifolds, then use this as a 
boundary condition for the solution on the 2-dimensional 
manifolds, which becomes in turn a boundary condition 
for the solution on 3-dimensional manifolds, etc. The 
leading-order rate requires only the solution for p = 0, 1. 
If we enumerate the lines connecting a point a/c G A to a 
point bk ^ B and use an arc-length parameterization for 
the kth line whose total length is then this given an- 
alytically as qo{s) = Q^^ J^^^ {ha{s^)I{s^))~^ ds\ where 



Qk= r {K{s')i{s')) ^ds' 



(28) 



is the normalization factor. On all other lines qo{s) = 0. 

Any dividing surface S hits each line at a single point, 
so the leading-order rate (in dimensional units) is 



<:ab 



AB 



(29) 



where the sum is over all connecting lines. This is asymp- 
totically equivalent to the rate one would obtain simply 
by restricting the full committor function and invariant 
measure to the set of 0- and 1-dimensional manifolds. 



To do this we took the set of 1-dimensional solutions 
computed as part of the free energy landscape, and com- 
puted the factor Qk from (29) on each manifold. Sum- 



ming over all of the 1-dimensional manifolds that connect 
a 0-dimensional manifold numbered a to a 0-dimensional 
mode numbered 6, gives the transition rate between a, 
h. We also include transitions between different ground 
states belonging to the same mode (e.g. a ^ a), by mul- 
tiplying the previous calculation by 2 since transitions 
can go in either direction along the line. 

Figure |6] shows shows the network of 0-dimensional 
states and the reaction rates between each state. The 
numbers reported are the dimensionless, purely geomet- 
rical parts of the rates Z^^ Qk^^ should be mul- 
tiplied by n~^D/d^ to give the dimensional rate. These 
rates, when multiplied by the total time of a simulation 
or experiment, give the average number of transitions one 
would expect to observe, so they are equal for both direc- 
tions a ^ b and b ^ a since our system is time-reversible. 
To obtain the rate relative to a particular state, i.e. the 
rate at which one leaves state a to visit state b next, given 
that the last state visited was state a, one should divide 
by the so-called reactive probability of a [53 . To leading 
order, this is equivalent to dividing by the equilibrium 
probability of a. 

Simulations We have verified our results by perform- 
ing Brownian dynamics simulations of ([2| with a short- 
range Morse potential. The results agree very well with 
our calculations of both free energy and transition rates. 
Figure [7| shows a comparison of the simulated proba- 
bilities versus theoretical probabilities of each mode for 
n = 6 (see SM for n = 7,8), for particles interacting 
with a Morse potential with range parameter p = 30, a 
range of ~ 5% of the particle diameter. The agreement is 
nearly perfect. Fig. [7] also compares the number of each 
type of transition we saw in the simulations, to that pre- 
dicted from theory. The theory slightly overpredicts the 
total number of transitions - this is what one should ex- 
pect from the geometrical picture, as these leading-order 
rates neglect the time the system spends in the floppy 
manifolds, which would tend to slow it down. 

These results are encouraging partly because they are 
evidence that we have executed these calculations cor- 
rectly, but also because they suggest the asymptotic limit 
may apply for experimental systems, such as Meng et. al. 
[4 where the potential reportedly had roughly this width. 



Comparison with other numerical approaches to the 
free energy landscape 



Transition Rates for hard spheres 

We used the formalism to compute the leading-order 
rates for n = 6,7,8 hard spheres with diameter d = 1. 



We have compared our results with those obtained by 
a numerical study that directly searched an energy land- 
scape of a short-ranged potential (a Morse potential, with 
range roughly 0.05 particle diameters) for local minima 
and transition states ^42^. The numerical method found 
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fewer local minima than there are 0-dimensional modes, 
and fewer transition states than 1-dimensional modes, 
suggesting that the asymptotic theory has the roughest 
landscape. By computing adjacency matrices for the nu- 
merical states with a cutoff bond distance of 1 + 1 x 10~^, 
we verify that for n = 6, 7, 8 each local minimum corre- 
sponds to a unique 0-dimensional mode, and each tran- 
sition state lies on a unique 1-dimensional manifold. 

We identify the point on the 1-dimensional manifold 
that is closest to each transition state. The transition 
states are very close to the local maxima of the vibra- 
tional factor —logh{x) (see Appendix for plot), consis- 
tent with it being a saddle point of the potential energy. 
We believe the small discrepancy in location can be at- 
tributed to the finite width of the Morse potential used 
in the numerical procedure. 

The missing 0-dimensional modes occur for manifolds 
that are very close to each other in the quotient space 
metric (separation ^ 0.08), which is the case for modes 
{1,4} (n = 7) and modes {1,2,3,4} (n = 8). For these 
only one local minimum was found for the entire group. 
We hypothesize that because the separation is within the 
range of the potential, the 0-dimensional modes merge to 
form a single local minimum. 

The missing 1-dimensional modes often, but not al- 
ways, correspond to self- self transitions - these transi- 
tions do not matter when particles are identical, however 
they will account for transitions between different states 
when the particles are not all the same (e.g. [54 ). For 
example, for n = 8 the numerical procedure identifies 
45 transition states, compared to our 75 one-dimensional 
manifolds. Of the missing manifolds, 16 are self-self tran- 
sitions, 9 are nonself-nonself transitions within the group 
{1,2,3,4}, and 5 are nonself-nonself for endpoints not 
both in the group. 

DISCUSSION / CONCLUSIONS 

We have developed a new framework for understanding 
energy landscapes when particles interact with a short- 
ranged potential. We show that in the limit as the 
range goes to zero and the depth goes to — oo, the en- 
ergy landscape becomes entirely governed by geometry, 
with a single parameter encapsulating details about 
the potential and temperature. When is large, only 
the lowest-dimensional geometrical manifolds contribute 
significantly to the landscape and this makes a computa- 
tional approach tractable. To illustrate the limit, we have 
computed the set of low-dimensional manifolds for n < 8 
hard, spherical particles. This solution to a nontrivial 
problem in statistical mechanics can be used to compute 
equilibrium or non-equilibrium quantities for any poten- 
tial whose range is short enough. 

We were able to calculate this set of low-dimensional 
manifolds because we began with the set of rigid clusters. 



and made the conjecture that all fioppy modes can be ac- 
cessed from these by breaking bond constraints. Solving 
for the complete set of rigid clusters is a difficult problem 
in discrete geometry that has only been done for n < 11 
|38l [55] , but with current computational power and novel 
approaches [55l |56] one can anticipate reaching larger n. 
Very large n will eventually require making approxima- 
tions to the geometry problem. We speculate that as n 
increases, structures with extra bonds, as well as "sin- 
gular" structures whose Jacobians have extra zero eigen- 
values, will come to dominate the landscape - these have 
not yet been considered in our asymptotic framework but 
they are observed with high probability in experiments 
[4. 

We have compared our results to those obtained by nu- 
merically searching the free energy landscape of a short- 
ranged Morse potential for local minima and transition 
states. Our method finds more minima and transition 
regions of the potential energy than the numerical search 
procedure, and this points to a potentially useful exten- 
sion of our theory - one can imagine starting with the lim- 
iting geometrical manifolds, and following these in some 
way as the range of the potential is increased, to obtain a 
low-dimensional approximation to the free energy land- 
scape for finite-width potentials, such as Lennard- Jones 
or Van der Waals clusters. This method would overcome 
a major issue with numerical searches which is that there 
is no way to ensure that all important parts of the land- 
scape has been found - we claim that our manifolds are 
the complete set of low-energy states so they will remain 
so under small enough perturbations. In addition, this 
would provide a way to deal with the increasing rugged- 
ness of energy landscapes with short-ranged potentials, 
which are a challenge for numerical methods - we start 
with the most rugged landscape and would only need to 
smooth it. 
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FIG. 1. Top: schematic of a traditional depiction of an en- 
ergy landscape. Local minima are separated by energy barri- 
ers (red line) that govern the transition rate from one basin 
to another. Middle: schematic of a geometrical energy land- 
scape, showing the 0- and 1-dimensional sets. Local minima 
are infinitesimally narrow, deep points, separated by long, 
nearly flat lines. Along these the dynamics is governed mainly 
by diffusion, so the length of the line (in red) determines 
the transition rate. Bottom: Example of a 1-dimensional 
manifold from the landscape for n = 6, showing the tran- 
sition path from polytetrahedron (left) to octahedron (right). 
Black line is the free energy Foc/ksT = — logriahal — 11 logK, 
along the 1-dimensional manifold in units of /c^T, where we 
chose K = 20. Black dots are the free energy Fa/ksT — 
— logUahal — 121ogAi: for the 0-dimensional endpoints. Red 
crosses mark locations of clusters that are plotted explicitly 
(they are plotted with 1/2 their actual diameter for clarity), 
and the horizontal axis is a parameterization of the manifold 
in the quotient-space distance s. 
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FIG. 2. Free energy landscape for n=6 with 0, 1, and 2 bonds broken. Red circles are 0- dimensional modes, blue are 1- 
dimensional, yellow are 2-dimensional. The area of each circle is proportional to the geometrical partition function riaCa of 
each mode, hence to its probability in equilibrium relative to modes of the same dimension. Modes are identified by numbers 
and arrows show the connectivity: an arrow from mode i to mode j indicates that mode i is part of the boundary of mode 
j. The number on each arrow indicates the number of different manifolds of type j that are connected to a single manifold of 
type i. The computed parameterizations are shown for each of the 2-dimensional modes. 




FIG. 3. A two-dimensional manifold (mode 18, n = 6), pa- 
rameterized in the plane, with selected points on the man- 
ifold plotted as clusters. The vertex of each triangle repre- 
sents a cluster and black or red dots indicate the ones that 
are plotted. The corners (black dots) are rigid structures, or 
0-dimensional manifolds. The edges are 1-dimensional man- 
ifolds and points on these are obtained from rigid structures 
by breaking one bond, while points in the interior form a 
2-dimensional manifold and are obtained by breaking two 
bonds. The 1-dimensional manifolds, beginning at the oc- 
tahedron (left) and moving clockwise, are 7,5,5,7. The text 
indicates the type of bond that breaks or forms as one moves 
along each edge. 
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temperature temperature temperature 



FIG. 4. Relative yields of 0-dimensional, 1-dimensional, 2-dimenisional modes (neglecting all higher- dimensional modes) 
for n — 6, 7, 8. The yield for p-dimensional modes is calculated as yp — t\j ^-^Zp/{h.^Zq + kZi + Z2) with K.{T) = 
g-(t/o/fcB)/T/^(^^//(^Q)/^^)(^2/7r)/T. We used [/o//cs = -4, ([/''(0)/A;B)(2/7r) = 15, but the numbers don't change the qualita- 
tive shape. Note that modes with dimensions > 2 will become important at the higher temperatures. 



P(x,y,t) 

i i 




p = p + p /x 



FIG. 5. Schematic showing the asymptotic boundary con- 
dition near a one-dimensional boundary. On the left, grey 
shading indicates the depth of the potential, and dashed line 
indicates the boundary at which the outer and inner solu- 
tions. By considering the total probability flux in/out of each 
small volume element near the boundary (red box), one can 
replace the detailed dynamics near the boundary in the limit 
as ere with an effective boundary condition (right). 
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FIG. 6. Rates (geometrical components) at leading-order, for n = 6,7,8. Dimensional rates are k, ^D/d^ times the above. 
These rates indicate the total number of each type of transition one expects to see, per unit time. 
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Theoretical counts 
Mode 1 2 

1 1570 ± 78 153 ± 24 

2 153 ±24 



Simulation counts 
Mode 1 2 

1 1256 124 

2 124 



1-d modes.Time = 2326 
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2-d modes, Time = 2326 
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FIG. 7. Figures: Theoretical probabilities of each mode 
(relative to modes of the same dimension), versus probabili- 
ties computed from simulations using a potential with width 
^ 5% of particle diameter, for n = 6. The red line is The- 
ory=Simulation indicating a perfect match. Tables: number 
of transitions of each type, for simulation and theoretical pre- 
diction. The theoretical prediction indicates expected 95% 
confidence intervals, computed using a normal approximation 
to a binomial. The sticky parameter was At = 16 and the total 
running time was 2.3 x 10^ units. 
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Parameterization of a neighbourhood of Qcx 

We provide a brief argument for the following state- 
ment made in the main text: given a regular point 
X G fta^ there exists a differentiable parameterization of 
a neighbourhood Af{x) C M^^ of the form y x {^a-l^i, 
with y e R^^-^, such that \/y • Vy«, = on 1^^.. 

Recall that a regular point x G l^c^ is a point 
such that the Jacobian of the transformation x 
{ya^{x)^ . . . ^ya^{x)) has rank m. If x is regular then 
it has a neighbourhood Ma{x) C Vta that is a differen- 
tiable manifold with co-dimension m ^\ ES], so there 
exists a parameterization near x hy y ^ ]^3n-m. ^qx^ ^.j^^ 
associated mapping be z : M^"^-^ Q,^, 

Given some point x' G (not necessarily on Vta)i 
we define a mapping x' {Va^.y) as {ya,{x'),y{x')), 
where y{x') is found from the limit of the gradient flow 
map, i.e. y{x') = ^""'^(limt^oo ^(^)) where (j) solves 
f = -"^EiUiVaM), m = x'. (We abbreviate 
to mean the full list of constraint variables.) Results 
from [59] (see also [60]) show that this map exists and is 
smooth enough in a neighbourhood C 1^ of x G l^a, pro- 
vided is regular and U {y) is sufficiently smooth. Since 
the mapping x' {Vai-, y) has full rank at x G ^1, it does 
also in a neighbourhood M{x) C and so by the Inverse 
Function Theorem it is invertible. The orthogonality at 
Via follows because this is a level set of the constraints. 

Note that while this provides the required parameteri- 
zation in a neighbourhood of x G I^q,, we have not shown 
that it extends to the set (see ([7|)). This is not a 
problem for the asymptotic calculations, as these remain 
valid if Vt^^ is replaced with an atlas of local parameteri- 
zations A/'(x), patched together with a partition of unity 
- the asymptotics only require the local behaviour near 
ycy. = and are not sensitive to the cutoffs at Tc. 

Quotient space metric and equations 

In this section we give more details about the quotient 
space and metric structure on it, and show how these 
arise naturally from our equations. Although these facts 
are well-known in Riemannian geometry [58l[61] and well- 
used in chemistry and mechanics [40l EH |63] , we have not 
found a reference dealing succinctly with our particular 
context so we collect the relevant facts and demonstra- 
tions here. 

The manifold structure on the quotient space Recall 
that we defined the quotient space associated with a man- 
ifold fla to be n2 = ^a/G, where the G = SE{3) is the 
Special Euclidean group, i.e. the group of rotations and 
translations of a cluster. Then is a smooth manifold 
if the Lie group G acts properly and freely on ([62l, 
Prop. 4.1.23 p. 266, [61]). To act properly is a compact- 
ness condition and it can be checked that it is satisfied 
for SE{3). To act freely means that the only element 



g G G such that g • x = x is the identity. This is the case 
provided there is no cluster in such that the spheres 
all lie on a line; for floppy manifolds with up to two bonds 
broken this is true when n > 5. 

The orbit of a cluster x is the set of points of the form 
g • X ioT g G G, and is written as [x]. Each orbit is 
identified as an element in Vt^ by the canonical projection 

Ti'.Vtoc^ ^a/G := {[x] : X G Qa}- (A.30) 

This projection shows how to map the tangent spaces 
to each other, via the pushforward map. Let 7^(x), 
T^{[x]) be the tangent spaces at x G 1^^, [x] G ft^ 
respectively. The tangent vectors map as follows: if 
c{t) G is a curve such that c(0) = x, then the tan- 
gent vector c'{0) G Ta{x) maps to the tangent vector 
^7r(c(t))|t=o ^ '7?([^])- Note that tangent vectors can 
also be identified as derivations, which we will write as 

Metric on the quotient space The group G acts iso- 
metrically on which means it respects the inner prod- 
uct (•, •)^^ on the manifold: given ^1,^2 G 7^(x), the in- 
ner product satisfies (^1,^2)^^ = {g ■ ti^g ■ ^2)^^ for all 
g G G. Therefore one can construct a metric on the 
quotient manifold that is compatible with the projection 
(making the projection a Riemannian submersion)^ and 
this metric is unique [6T] . 

To specify the metric we decompose Ta{x) into the ver- 
tical subspace T^{x), containing the directions tangent to 
the action of G, and the horizontal subspace 7j(x), its or- 
thogonal complement. Therefore 7^ (x) = {x)^Tf^{x). 

Let Pa{x) : Ta the the orthogonal projection 

operator (we omit the argument x for succinctness.) Let 
[t] denote an element of that has representative t G 
Ta- The metric g^ on the quotient space Vt^ is computed 
from the metric g^ on Q.^ as 

([^1], M)gQ = {VatuVat2)g^- (A.31) 

Fokker-Planck equations on the quotient space Next 
we show how this differential structure arises naturally 
as a result of our manipulations to the Fokker-Planck 
equation. 

Consider a point xq G fta^ and let us parameterize a 
neighbourhood on ft^ in such a way that at Xq, the direc- 
tions tangent to infinitesimal rotations and translations 
are orthogonal to the remaining variables. This can be 
done with the standard Euler angles. Let 

/I \ 

Rx{0) = COSI9 - sin i9 
\0 sin i9 cos i9 / 

(cosO sin^\ 
10 
- sin i9 COSI9/ 

/ cos 6 — sin 0\ 
R^{0) = sinO cosO 
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be the matrices for rotation of a point about the x^y^z 
axes respectively, with block-diagonal versions appropri- 
ate to a cluster of n particles Ra;(6>), ^z{'d)- These 
are obtained as, for example. 



\ 



with n copies along the diagonal and zeros everywhere 
else, and similarly for the other matrices. Let 

T(/ii,/i2,M3) = (/il,/i2,M3,Ml,-.-,/^l,M2,M3)^ 

be a vector representing translations. Let (p^y) : MP 
R^^ parameterize the remaining directions with ^(0) = 
xo, so that a neighbourhood of xq on fta can be param- 
eterized with variables {Ox^ ^y, Mi, M2, Ms, y) as 

X = Il,{0,)Ily{Oy)Ilx{Ox)Hy) + T(/ii,/i2,/i3). (A.32) 

We can choose (j){y) so that it is orthogonal to infinitesi- 
mal rotations and translations at xq: this means that we 
\y=o ' L,=o^o) = ' each yj, and each 



require 



dyk I 



rotation matrix; the condition for translations is satisfied 
if the center of mass of xq is at the origin. 

The metric tensor on is Qa = J where J is the Ja- 
cobian of the transformation (^a^, 6^^, Mi, /i2, Ms? ^) ~^ 
X. This has a simple block diagonal structure at xq: 




(A.33) 



where 



XjZi 



-yiZi . (A.34) 



Here ga is the pxp contribution from the ^-variables. Is 
is the 3x3 identity matrix, and I is the moment of inertia 
tensor, which assumes the configuration is written as x = 
(xi, yi^ Zi^ . . . ^Xn^yn^ ^n)- We will write its determinant 
as P{xo) = det(I). 

Therefore the Fokker-Planck equation on after in- 
tegrating over the fast variables has the following form 
at Xq: 



13a 



(A.35) 



Here o,,-, 9uv are the elements of ^q, , ^q, respectively, 
a, h index the rotational and translational variables, and 



we have substituted (A.33) in the second equation. We 



I 

that we chose for [xq] because G acts isometrically. 
After identifying functions, tangent vectors, and the 



will not deal with the fiux explicitly as this comes from 
the same manipulations on the higher-dimensional man- 
ifolds. 



Integrating (A.35) over orbits gets rid of the second 



term on the RHS, by Stokes' theorem, so we are left with 
1 



dt{K^Kip) 



'Kir^'d.p)- 



(A.36) 



where p is the integrated value of p on an orbit (this is Cp 
if p is constant on the orbit, where the constant C does 
not depend on xq.) Because this does not depend on the 
location along the orbit of a point, we can identify it with 
a function on the quotient space as p{x,t) = p^{[x],t). 
For the same reason we can identify tangent vectors via 
the canonical projection as duP = 9[n]P^- The metric 
guv has the same elements as the quotient metric g[u][v] 
because it only involves tangent vectors in the horizontal 
subspace, and it is independent of the representative xq 



metric in (A.36) with their projections in the quotient 
space, we obtain (20) in the text. 

Representing the quotient space To parameterize the 
quotient manifold it is convenient to map it to a space 
that has an explicit representation. For our numerical im- 
plementation we store the edge-lengths of all the particles 
(we call this "bond-distance" space.) This representation 
* actually forms the quotient space with reflections as well, 
so it is only diffeomorphic to does not contain 

a cluster where all of the particles lie in a plane. 

An alternate parameterization would be to constrain 
one vertex to be at the origin, one vertex to lie on the 
X-axis, and one vertex to lie on the x7/-plane. This would 
embed in Qa- 

Given a parameterization of the quotient manifold in 
some space B with canonical projection tt, the tangent 
vectors and metric can be computed from the pushfor- 



ward map (A.30) via numerical differentiation. That is. 



given a point [x] G B with representative x ^ fta where 
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this is two-dimensional, we compute the two unit tangent 
vectors ^1,^2 ^ Ta{x) that are perpendicular to the in- 
finitesimal rotations and translations. This is easy to do 
from the null space of the matrix M defined in section [] 
We take small steps in each of these directions to obtain 
points xi = X + As ti, X2 = X -\- As ^2, project to and 
obtain first-order estimates of the quotient tangent di- 
rections as [ti] = (7r(xi) —x)/As, [^2] = (7r(x2) —x)/As. 
These have lengths 1 and inner product (^1,^2)^^, which 
defines the metric on B. 

We use a first-order scheme to find the distance be- 
tween two nearby points [xi], [X2] in B. We find the sep- 
aration vector V = [xi] — [^2], project this onto the tan- 
gent space at [xi] and find the length of this projection. 
We repeat at the tangent space to [^2] and average the 
lengths. 

Note that parameterizing the quotient manifold as a 
subset of would imply slightly different numerical 
algorithms. For example, to find the distance between 
two nearby points, one would simply project the separa- 
tion vector V onto the horizontal tangent space at xi, X2 
(these representatives can be chosen equal to [xi], [X2]) - 
this might be more efficient than the bond-distance space 
method. 



Adjoint equations 

In this section we compute the backward Fokker- 
Planck equation associated with ( [2Q| ). We could derive 
this using the same asymptotic procedure on the back- 
ward equation associated with ([2|, however we prefer to 
demonstrate how to convert between the forward and 
backward sticky equations directly. We only outline the 
arguments here, leaving several steps to the reader. 

Recall that the backward equation describes the evolu- 
tion of the expected value of some function g{x) 
that starts with a unit mass at x and is subsequently 



stirred by the probability dynamics. [6T This is obtained 
from the transition probability measure Px{dy^t) as 



Px{dy,t)g{y) = / p{x,y,t)g{y)dp{y) 



u{x^ t) - 

(A.37) 

where the transition probability has initial condition 
Px{y^O) = "^^Saiy — x) and we have used the fact 
that it has a density y^ t) with respect to the equi- 
librium measure dp{y) (see ([24|)). Note this implies 

p{x, y, 0) = Y.a ^cx{y - X)/K.a' 

We suppose that Px satisfies the Chapman- 
Kolmogorov equation as this property holds for the 
original probability measure p^dx from which it derives: 



L 



P^{dz,t - s)P,{dy,s) = P^{dy,t). 
This allows us to write 



(A.38) 



u{x 



^t) = Px{dz,t-s)u{z,s) = / px{z,t-s)u{z, 



s)dp{z). 



(A.39) 

We can now obtain an evolution equation for u. Applying 



dt to ( |A.39[ ) and using ( |24[ ) gives 
PtU dp{y) = 



du 

m 



( div grad p)u dp{y) 



= p{ div grad u) + / gi"ad p ■ n) 

— p{ grad u ■ n) -\- hZiU div grad p. 

Here {fli} is the set of manifolds of co-dimension 1 that 
form the boundary of the full space 1], hhri are the sticky 
factors along these manifolds, div , grad denote dif- 
ferential operators on Q and div^ , grad^ will denote 
those on h denotes a generic outward normal to the 
appropriate manifold, and integration is with respect to 
the volume element appropriate for each manifold. We 
now substitute for grad p • n using the forward sticky 
equations to obtain 



du 



/ Pi div grad + V / «( div, grad, p - d.v grad p) - p grad u ■ n + ..u d.v grad p 
I Pi div grad u) + y:[ Pi div, grad, . - grad . • n) + V / V grad, p ■ n(^^) - p grad, . • n('^)) 



(A.40) 



Here {^j} is the set of manifolds of co-dimension 2, form- 
ing the boundaries of the manifolds ; the sum in the fi- 
nal term is over all manifolds Qi that have as a bound- 
ary, and n^*-^^ is the outward normal vector from fti at 



Evaluating (A.40) ai s = t gives the backward sticky 
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equations 

dtu = div grad u in Q 

hCidfU = diYi Ki grad^ u — Vu • n inVti (A. 41) 

dtU = (boundary terms) in Vtj 

We stop at the first two terms; the interested reader can 
show that evaluating the boundary terms leads to the 
expected equations on the lower-dimensional manifolds. 
This set of equations is identical to the forward equations, 
so the system is self- adjoint. 



Parameterizing the manifolds 

In this section we outline the method we used to pa- 
rameterize the quotient manifolds O!^ describing clusters 
of hard spheres with up to 2 bonds broken. The method 
can be broken down into two separate sets of algorithms. 
The first algorithm generates points on the manifold, by 
taking linear steps along the tangent directions and pro- 
jecting back down to the manifold. This is sufficient to 
calculate the 1-dimensional manifolds. The second algo- 
rithm links up the points with bars to form a simplex on 
which calculations can be performed, and is required for 
2- and higher-dimensional manifolds. 



1-dimensional manifolds 

Consider first the 1-dimensional manifolds. We take 
steps along the manifold as follows: given a point xq G 
Qq,, a set of bond-distance constraints {yk}^^i-, and a ba- 
sis {ti}^^i for the part of the tangent space at xq parallel 
to rotational and translational motions, we form a ma- 
trix M = (V^i, . . . , V^rn, ^1, • • • , te)^ and compute the 
null space of M. When xq is a regular point on the man- 
ifold, this null space contains a single vector v lying in 
the tangent space of Qa{xo)^ so we take a step in that 
direction as xi = + {As)v. Because v is orthogonal to 
translations and rotations, the length of our step in the 
quotient metric is | |xi — xo| | J = As + O(As^). 

This step pushes us slightly off the manifold, so we 
project back down to it by finding a set of Lagrange 
multipliers so that the projected point x[ = Xi -\- 
"^Zk ^k^yki^i) lies on the manifold, i.e. we solve the non- 
linear system of equations yk{xi + ^k "^k^Vki^i)) — 
|65 . This is easily done using Newton's method as we 
are typically very close to the manifold. 

Beginning with a rigid cluster xq with an associated 
set of bond constraints, we break a bond by deleting one 
of the constraints, and perform the steps above until an- 
other bond is formed. This provides an ordered set of 
points in the quotient manifold, along with distances be- 
tween them - this is a "line" . 



2-dimensional manifolds 

Computing the 2-dimensional manifolds is more in- 
volved; here are the steps we followed. 

Compute the boundaries First, we compute the 
boundaries of the manifold. To compute a 2-dimensional 
manifold with constraint list yi{x), . . . ,yn{x), we 
start with a corner point Xq G ^2 with two extra con- 
straints yi^{x),yi^{x). Deleting one of these, say yi^, we 
walk along the 1-dimensional manifold (as in the previ- 
ous section) until another bond is formed, corresponding 
say to constraint yi^. We add this to our constraint list, 
delete the next extraneous constraint yi^^ and repeat. We 
continue in this way, deleting one extraneous constraint 
at each corner, until we reach the original corner xq. This 
gives us the boundary of including corners. 

Generate points in the interior Second, we generate a 
collection of points in the interior. Beginning with every 
point on the boundary, we generate lines in the interior by 
holding all constraints fixed except the one that takes us 
off the boundary, and walk in this direction until we exit 
the manifold. We throw away points that are too close 
in some metric to existing points to avoid generating too 
many points. 

Triangulate the points Third, we triangulate the 
points. Many typical algorithms will not work here be- 
cause our surface is not embedded in R^, so we adopt an 
algorithm proposed by by [66] (see also [67 ) that works 
as follows. 

1. Map the boundary to a fixed convex polygon in M?. 

2. Map the interior points to the interior of that re- 
gion in R^, by letting each interior point be a con- 
vex combination of its neighbouring points. More 
specifically: let xi^...Xn be the set of interior 
points, and let Ni be a neighbourhood of an interior 
point Xi. Given a set of strictly positive weights A^j 
such that Xla^jGA^i ~ -^^^ parameter points 
ui^ . . .Un G R^ that solve the linear system of equa- 
tions 

Ui= XijUj, i = 1, . . . ,n. 

XjeNi 

Each Ui is contained in the convex hull of its neigh- 
bours, so it will be in the interior of the region 
defined in step (1) [66 . 

3. Triangulate the parameter points in R^ (we use a 
Delauney triangulation.) This lifts back to a trian- 
gulation of the manifold. 

To implement this, we choose the boundary region so 
that the corners lie on a circle with a fixed radius, and 
the line segments joining them lie on arcs of circles whose 
lengths are approximately the same as the lines they are 
parameterizing. Points along these arcs are placed so the 



22 



inter-point distance in the plane is proportional to the 
distance in the quotient metric between the points. The 
quality of the triangulation will depend on the choice 
of boundary region, and we find better qualities as the 
angles at the corners more closely represent the angles 
on the manifold. 

Choosing the neighbourhood Ni is a balance between 
sampling many points to get a smoother parameteriza- 
tion, and choosing fewer points so the manifold is roughly 
linear in the neighbourhood and does not contain any 
folds or other external branches of the manifold. We 
choose the neighbourhood to be the k nearest points 
along the manifold, where a range of roughly 8 < k < 15 
works well for the step sizes we use, but k will increase 
as step size decreases. 

There are many ways to choose the weights; the most 
straightforward is for them to be the same, but almost as 
straightforward is for them to be inversely proportional to 
distance in some metric. For rapid but still good quality 
results we use the metric in bond-distance space - this 
takes the list of pairwise bond distances and computes 
the Euclidean distance between the vectors. 

Improve the triangulation Finally, we improve the 
quality of the triangulation by letting the triangle sides 
be springs and evolving the points on the manifold 
with the spring forces, re-triangulating when necessary. 
Springs are chosen to be slightly longer than the average 
distance between points so the points want to spread out 
and fill the whole space. When a point hits a boundary 
it is absorbed, and is subsequently constrained to move 
along the boundary. This algorithm was introduced by 
[68^ for a triangulation of the Euclidean plane, and we 
adapt it by replacing the length in the plane with the 
distance metric in the quotient manifold. In practice, we 
typically use bond-space distance instead as this is faster 
to compute and still gives a good quality triangulation. 

Integration on the manifolds To compute quantities 
integrated over the manifolds, such as for the partition 
functions z^^^^^ we used finite elements on the simplex 
with standard piecewise linear elements. The sides of the 
triangles must be calculated in the quotient space metric. 

Remarks 

Topology The calculations above require that the 2- 
dimensional manifolds be topologically equivalent to a 
disc. Because we have obtained smooth parameteriza- 
tions, we are confident that this is true for all the fioppy 
manifolds under consideration. Alternatively, one could 
show this from a collection of points using Betti numbers 
[69] , for example. 

This points to an interesting question in discrete ge- 
ometry - under what conditions is the topology of fioppy 
manifolds a polygon? One has to rule out surfaces of 
higher genus and non-orient able surfaces, among other 



things. We expect this can be shown for low-dimensional 
manifolds and small n, but larger p or n may be more 
complicated. 

Numerical parameters and convergence For the cal- 
culations reported in the text we used a step size of 
As = 0.01 to calculate the one-dimensional manifolds, 
and As = 0.05 to generate points on the boundary and 
interior of the 2-dimensional manifolds. When generating 
points we threw away points that were closer than 0.5 As 
in bond-distance space to already- generated points. We 
typically ran the triangulation step 3 times before re- 
triangulating, using a step size of At = 0.1 and an in- 
ternal pressure parameter of 1.2 (see [68 ), although a 
small selection of manifolds had to be re-triangulated 
after each spring step. We ran the triangulation un- 
til the area of the manifold, computed in bond-distance 
space, changed less than A5/2O after 3 consecutive re- 
triangulations. We checked for convergence in two ways: 
by calculating the manifolds using a coarser resolution 
{As = 0.1 for the 2-dimensional manifolds. As = 0.05 
for the 1-dimensional manifolds), and by running the tri- 
angulation algorithm for longer, until the minimum tri- 
angle quality {q = 2rin/routi the ratio between (twice) 
the radius of the largest inscribed circle and the small- 
est circumscribed circle, see e.g. [70 ) was greater than 
0.2. Both tests allowed us to conclude that our calcu- 
lated ratios Z2/Z1, Zi^Zq are correct to ±0.05, ±0.01 
respectively (although we have not reported this many 
decimal points for the latter). 



Simulations 

We have performed Brownian dynamics simulations 
of interacting particles to test our asymptotic calcula- 
tions. We solve equation Q with I) = /3 = 1 us- 
ing a forward Euler timestep. For the potential we 
started with a Morse potential with maximum depth E 
at r = 1 and range parameter p; this takes the form 
U{r) = Ee-P^"^-^^ (e-^^^-^) - 2). The hard-core part for 
r < 1 was modelled with a parabolic potential of the 
form U{r) = ^m?U" {l){r — 1)'^ — E for some number m; 
we choose m = 2. This constrains the time step to be 
At <C {m?U" and we typically choose a factor of 
6-8 less. The sticky factor, accounting for the parabolic 
part, is = ^^^^A/f • 

The whole potential was truncated at Tc = 1 ±4/p, by 
adding a linear term to keep the force continuous at Tc, 
as Utruncir) = U{r)- {U{rc) ^U\rc){r - rc)) for r < r^ 
and Utrunci'^) — otherwise. This modifies the sticky 
parameter to 

m ± 1 exp {E - Ujrc) - U'{rc){l - r^)) 
^ m V2Ep^ V 2* 

(A.42) 
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We began with an initial condition drawn from the 
equihbrium distribution of rigid modes, and ran several 
copies of the simulation for a very long time. At time 
increments of 1 x 10 ~^ we checked to see which floppy 
mode the cluster occupied by forming its adjacency ma- 
trix, computing the number of bonds, and if this number 
was > 3n — 6 — 2, finding the adjacency matrix in our list 
of fioppy modes to which it was topologically equivalent. 
This allows us to compute the occupation probabilities 
of each fioppy mode. 

To form the adjacency matrix, we said particles were 
bonded when they were at a distance of less than 1 + 2/p; 
this is the range beyond which the force is negligi- 
ble. (Choosing a bond distance of 1 + 1/p gave results 
that were inconsistent with the asymptotics.) We used 
the Matlab function graphisomorphism.m to determine 
topological equivalencies. 

We also kept track of transitions between rigid clus- 
ters, by recording the times and mode numbers at which 
the system first hits a rigid state that occupies a sepa- 
rate part of configuration space than the previous rigid 
state. This is easily done by checking whether or not the 
current adjacency matrix of a rigid mode is identical to 
the adjacency matrix of the previous rigid mode. 

Figure [S] shows the simulated probabilities versus the- 
oretically computed probabilities of each fioppy mode for 
n = 7, 8, for simulations using a Morse potential with 
the same parameters as in the text (Figure [t]): E = 8.5, 
p = 30, dt = 2x 10~^. This implies the sticky parameter 
is K: = 16. Again, there is excellent agreement. 

Figure [9] plots the elements of the transition count ma- 
trix for n = 6, 7, 8 from simulations, versus two theo- 
retical calculations. The blue markers come from the 
leading-order asymptotic approximation, computed from 



O-d modes, Time = 6184 



1-d modes.Time = 6184 



(29). The magenta markers come simply restricting the 



dynamics to the network of points and lines, which gives 
geometric rates of (Zq + ^kQk^ ~ "these are 

the blue rates multiplied by -\- Zi/Zq)^ so are 
asymptotically equivalent, but uniformly smaller. These 
rates do a better job of predicting the simulated rates for 
this value of the range parameter. 

To calculate rates for n = 7, 8 we have grouped the 0- 
dimensional modes that are very close together into one 
mode (these are modes {1,4} for n = 7, modes {1, 2, 3, 4} 
for n = 8) . The rate from a mode in this group to a mode 
not in the group is the sum of the rates out of each mode 
in a group, and the rates within a group are ignored. We 
were able to make a distinction between modes in a group 
only by using a range parameter of p = 150 (not shown). 

Table [ll| shows the ratios Zp+i/Zp extracted from the 
simulations. These are uniformly smaller than the the- 
oretically computed values, but approach the theory as 
the range of the potential decreases. This is because the 
simulations use an excess bond length that is finite rather 
than infinitesimal, so some clusters are counted as being 
on p-dimensional manifolds, when in the limit they would 
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FIG. 8. Simulated versus theoretical probabilities of floppy 
modes, for n = 7 (top), n = 8 (bottom), for a Morse potential 
with range parameter p = 30. The running time of each 
simulation varied and is given in the title. 
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FIG. 9. Elements of the count matrix for n — 6, 7, 8 (left, 
middle, right). 



be on p+l-dimensional manifolds. One could potentially 
correct for this if one had no knowledge of the potential, 
but wanted to use these measurements to estimate the 
sticky parameter. 



Data 



We provide details for the fioppy modes for n 
Details for n = 7, 8 are available upon request. 



6. 



Table III reports the following quantities for each 



mode: the volume S = J^q 1 , mean geometrical sticky 
parameter h = {J^q h)/S^ mean rotational contribution 



C 



Zi/Zo 
A B 

6 5.5 5.9 

7 5.8 6.0 7.2 

8 5.2 6.4 6.7 



C 



Z2/Z0 
A B 

6 3.4 4.1 

7 3.5 3.8 4.2 

8 3.6 3.6 4.7 



TABLE II. Ratios extracted from simulations with different 
range and sticky parameters. These were computed as (time 
in p + 1-dimensional modes) / (t ime in p-dimensional modes) 
/ K,, where k, is computed using ( |A.42 ). Parameters were (A) 
dt = 2e - 6, E = 8.5, p = 30 V= 16), (B) dt = Se - 7, 
^ = 10, p = 50 (At = 31), (C) dt = 2.5e - 7, ^ = 10.6, 
p 150 {k 22) 
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FIG. 10. Vibrational contribution to the free energy —\ogh{x) along selected 
indicate the transition state computed in [42 for a potential with finite width. 



1-dimensional manifolds (red line). Markers 



Mode h 



Za corners 



0-dimensional modes 



1 


0.061 3.16 


180 


34.64 












2 


0.034 2.83 


15 


1.44 












1-dimensional modes 
















3 


0.066 3.30 0.85 


180 


33.30 


1, 


1 








4 


0.063 3.29 0.89 


90 


16.65 


1, 


1 








5 


0.057 3.29 0.95 


360 


64.03 


1, 


1 








6 


0.069 3.49 1.47 360 


126.89 


1, 


1 








7 


0.045 3.04 0.63 


180 


15.40 


1, 


2 








2-dimensional modes 
















16 


0.075 3.37 0.35 


180 


15.99 


1, 


1, 


1 






17 


0.075 3.76 2.00 


360 


202.45 


1, 


1, 


1, 


1, 1 


18 


0.083 4.01 2.17 


180 


130.10 


1, 


1, 


1, 


1 




19 


0.064 3.53 1.07 


360 


87.44 


1, 


1, 


1, 


1 




20 


0.057 3.17 0.23 


180 


7.52 


1, 


1, 


2 






21 


0.073 3.79 2.84 


360 


284.23 


1, 


1, 


1, 


1, 


1, 1 


22 


0.055 3.17 0.24 


90 


3.76 


1, 


1, 


2 






23 


0.064 3.56 1.55 


72 


25.33 


1, 


1, 


1, 


1, 


1 


24 


0.067 3.48 0.64 


360 


53.23 


1, 


1, 


1 






25 


0.063 3.59 1.58 


360 


129.53 


1, 


1, 


1, 2, 


1 


26 


0.054 3.27 0.59 


360 


37.56 


1, 


1, 


2, 


1 




27 


0.081 4.07 2.67 


120 


105.52 


1, 


1, 


1 






28 


0.072 3.77 2.39 


90 


57.97 


1, 


1, 


1, 


1 





/ = (J^Q I)/S, multiplicity Ua (divided by a constant). 
Each of the integrals is over a single manifold in the set of 
isomorphic manifolds, and the average for rigid modes is 
simply the point value. For floppy modes we also report 
the corners, in the order they occur as one travels around 
the boundary. 



Figure [TT] provides a way to identify the individual 
modes. It shows the rigid clusters with particles num- 
bered. The table indicates how to reach each floppy mode 
by starting with a given rigid structure and breaking se- 
lected bonds. In most cases there are multiple ways to 
reach each floppy structure and every possible bond com- 
bination that does so is listed. 



TABLE III. Free energy data for each mode, n = 6. 



Figure 10 plots the transition states computed in |42] 
along each of the identified 1-dimensional floppy mani- 
folds. 



25 



FIG. 11. Rigid clusters for n = 6. Right, Mode 1 (polytetrahedron), Left, Mode 2 (octahedron). The table below shows which 
bonds to break to reach each floppy mode. 





Mode Starting Cluster break bonds 



10 
11 

12 



13 

14 

15 
16 
17 



18 



19 
20 



1-3, 2-4 
1-4 

1- 5, 4-5, 

2- 5, 3-5, 
5-6 

1-3, 2-3, 



■6, 4-6 
■6, 3-6 



{1-3 
{1-3 
3-6} 
{1-3 
{1-3 
2-6} 
{1-3 
{1-3 

2- 6} 
{1-4 

3- 6} 
{1-4 
{1-3 
{1-5 
{1-5 
{1-5 
5-6} 
{1-3 
2-6} 
{2-4 

4- 5} 
{1-5 
{1-3 
1-6} 
{2-4 
4-6} 
{2-5 
{2-5 



1-4} 

1- 5} 
{2-4 

3- 5} 

4- 5} 
{1-6 

5- 6} 

2- 3} 
{3-5 

1- 5} 
{1-4 
5-6} 

2- 4} 
4-5} 

1- 6} 

2- 6} 
{2-6 

1- 5} 
{2-3 
4-5} 
{1-6 
4-6} 

2- 5} 
{2-3 

3- 5} 
{3-5 

2- 6} 

3- 6} 



1-4, 2-4, 1-5, 2-5, 3-5, 4-5, 1-6, 2-6, 3-6, 4-6 
2-4}, {1-3, 2-4} 

2- 5}, {1-3, 1-6}, {1-3, 2-6}, {2-4, 3-5}, {2-4, 4-5}, {2-4, 
{2-5, 3-5}, {2-5, 3-5} 

3- 6}, {2-4, 2-5}, {2-4, 2-5} 

4- 6}, {2-4, 1-5}, {2-4, 1-6}, {1-5, 2-5}, {3-5, 4-5}, {1-6, 



{1-3, 
{1-3, 
4-6}, 
{1-3, 
{1-3, 
2-6} 
{1-3, 
{1-3, 
4-5}, 



5-6} 
1-4}, 
{3-5, 



{1-4, 2-5}^ 
4-6}, {1-5, 

{1-4, 5-6} 



1-4}, 
4-5} 
3-6}, 
1-6}, 



{2-3 
{1-5 
{1-5 
{2-5, 
5-6} 

{1-3, 3-5}, 
3-6}, {1-4, 
{2-4, 2-6}, 

3- 6}, {1-6 
{1-5, 5-6} 
{1-3, 4-5}, 

4- 6}, {1-4, 



{2-4, 
1-6}, 
{2-5, 
{2-5, 



1- 6}, 
{3-5 

2- 6} 

3- 6} 



{2-3, 2-4}, {1-4, 2-4}, {1-5, 2-5}, {1-5, 1-6}, {2-5, 

3-6}, {4-5, 4-6}, {1-6, 2-6}, {1-6, 2-6} 

{1-4, 3-5}, {1-4, 4-5}, {1-4, 1-6}, {1-4, 2-6}, {1-4, 

3- 5}, {2-5, 4-5}, {1-6, 3-6}, {1-6, 3-6} 

{1-5, 2-6}, {2-5, 1-6}, {3-5, 4-6}, {3-5, 4-6} 

{2-5, 4-6}, {3-5, 1-6}, {4-5, 2-6}, {4-5, 2-6} 
{2-5, 5-6}, {3-5, 4-6}, {3-5, 5-6}, {4-5, 3-6}, {2-6, 

{1-3, 1-6}, {1-3, 3-6}, {2-3, 2-5}, {2-3, 3-5}, {2-3, 

1- 5}, {1-4, 4-5}, {1-4, 1-6}, {1-4, 4-6}, {2-4, 2-5}, 
{2-4, 4-6}, {1-5, 3-5}, {1-5, 4-5}, {2-5, 3-5}, {2-5, 

4- 6}, {2-6, 3-6}, {2-6, 3-6} 

{4-5, 1-6}, {4-5, 5-6}, {1-6, 5-6}, {1-6, 5-6} 
{1-3, 2-6}, {1-3, 4-6}, {2-3, 1-5}, {2-3, 4-5}, {2-3, 

2- 5}, {1-4, 3-5}, {1-4, 2-6}, {1-4, 3-6}, {2-4, 1-5}, 
{2-4, 3-6}, {1-5, 3-6}, {1-5, 4-6}, {2-5, 3-6}, {2-5, 
2-6}, {4-5, 1-6}, {4-5, 1-6} 



